function [model, accval]=variational_EM(data, MAXCOUNT, MAXESTEPITER, MAXMSTEPITER, MaxFun)


model=[];

LOWLIMIT=0.01;
convergence=1000000;
model=init_params(data);
countVEM=1;
model.k=data.k;
model.k_hat=data.k_hat;
maxvalue = -1000000000000000;

[value1] = cal_likelihood(model, data);

while (convergence> LOWLIMIT && countVEM<MAXCOUNT)
    
%     [value1] = cal_likelihood(model, data);
%     
%     
%     if (compareval(value1, maxvalue))
%         maxvalue = value1;
%         %disp('correct');
%     else
%         disp('Incorrect');
%         keyboard;
%     end
    
    
    model    = E_step(model,data, MAXESTEPITER, MaxFun, maxvalue, countVEM);
%     [value2] = cal_likelihood(model, data);
%     
%     
%     if (compareval(value2, maxvalue))
%         maxvalue = value2;
%         %disp('correct');
%     else
%         disp('Incorrect');
%         keyboard;
%     end
    
    
    model    = M_step(model,data, MAXMSTEPITER, MaxFun, maxvalue, countVEM);
    [value3] = cal_likelihood(model, data);
    
%     if (compareval(value3, maxvalue))
%         maxvalue = value3;
%         %disp('correct');
%     else
%         disp('Incorrect');
%         keyboard;
%     end
    
    
    convergence=100*abs((value3-value1)/value1);
    value1=value3;
    value(countVEM)=value1;
    disp('count from V-EM');
    countVEM = countVEM+1
    
end

figure(2), plot([1:countVEM-1],value,'b.-');


probassgn = exp(model.alpha_1)./repmat(sum(exp(model.alpha_1),2),1,model.k);
[indgarbage,ind]   = max(probassgn');

accval = cal_accuracy(model, data);

end
